Prostorni podaci u R-u - vežbe

Uvod

---

Prostorni podaci su namenjeni skladištenju, manipulaciji i analizi informacija o nekim prostornim fenomenima. Na primer, izmerena temperatura u jednoj tački može da predstavlja jedan prostorni podatak. Prostorni podatak je najčešće reprezentacija jednog prostornog fenomena u toj tački. U ovom slučaju to je temperatura, koja je prostorno-vremenski promenljiva. Prema tome, prostorni podaci opisuju prostorne fenomene i kao takvi, prilagođeni su prirodi i karakteristikama prostornih fenomena:

• Prostorni fenomeni mogu biti diskretni i kontinualni,
• Pod diskretnim prostornim fenomenima se smatraju fenomeni koji se mogu opisati primarnim geometrijskim primitivima,
• Diskretni prostorni fenomeni se obicno opisuju vektorskim podacima (lokacija, kuca, reka…),
• Kontinualni prostorni fenomeni su obicno podložni prostornim i vremenskim promenama i obično se opisuju rasterskim podacima.
• Kad kažemo vektorski podaci, obično mislimo na tačku liniju i poligon.

sf paket

---

Razvoj novih tehnologija i softvera otvorenog koda za rad sa prostornim podacima, kao i standarda u ovoj oblasti, zahteva jednostavniji i efikasniji pristup nad podacima. R paket sp je razvijen da proširi R funkcionalnosti klasama i metodama namenjenim prostornim podacima. Postoje opravdani razlozi za prelazak na „najsavremeniju“ sf klasu, koja omogućava interoperabilnost i efikasnije rukovanje prostornim podacima. Paket sf je R paket za čitanje, čuvanje, rukovanje i manipulaiciju prostornim podacima u R okruženju, putem pristupa koji se naziva simple feature. Simple feature je međunarodni standard za predstavljanje i kodiranje prostornih podataka, predstavljenih putem geometrije tipa tačka, linija i poligon (ISO, 2004). Najvažniji razlog za upotrebu ovog pristupa je bolja manipulacija geometrijom prostornih podataka i jednostavan način za čuvanje geometrije u podacima, tj. jednom redu data.frame-a , čime se u potpunosti mogu iskoristiti funkcionalnosti familije paketa tidyverse.

Zadatak - Prostorizacija vrednosti emisije zagađenja na području Republike Srbije

---

Tekst zadatka

Vrednosti emisije zagađenja po polutantima SO2, NOx, PM10, PM2.5, NMVOC i NH3 prostorizovati po izvorima zagađenja (tačkastih, linijskih i površinskih kategorija). Vrednosti sumirati po ćelijama grida za područje teritorije Republike Srbije.

• Učitati u R okruženje izvore zagađivača sa vrednostima po svim polutantima. Podaci su u vektorskom obliku, u formatu ESRI shapefile.
• Učitati u R okruženje jedinstveni grid za teritoriju Republike Srbije, prostorne rezolucije 0.05° x 0.05°.
• Uraditi presek geometrijskih entiteta sa gridom.
• Sumirati vrednosti polutanata po ćelijama grida.
• Kreirati jedinstveni grid u kome ćelije grida imaju agregirane vrednosti po svim kategorijama.
• Uraditi vizuelizaciju podataka.
• Grid sačuvati u geopackage formatu za razmenu prostornih podataka, kao i gridovane podatke po kategorijama kao posebne lejere.

Podaci su dobijeni za sledeće kategorije:

• 1A1a - Public electricity production (tačkasti izvor)
• 1A3bi_R - Road Transport: Passengers cars - Rural transport (linijski izvor)
• 3Da1 - Inorganic N-fertilizers (površinski izvor)

Učitavanje podataka

---

Učitavanje tačkastih izvora zagađivača i transformacija:

sf.1A1a_ep <- st_read(dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/1A1a_ep.shp")

## Reading layer `1A1a_ep' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\1A1a_ep.shp' using driver `ESRI Shapefile'
## Simple feature collection with 21 features and 14 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 353634.2 ymin: 4795824 xmax: 587823.5 ymax: 5105985
## projected CRS:  WGS 84 / UTM zone 34N

sf.1A1a_ep 

## Simple feature collection with 21 features and 14 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 353634.2 ymin: 4795824 xmax: 587823.5 ymax: 5105985
## projected CRS:  WGS 84 / UTM zone 34N
## First 10 features:
##         NOx        SO2        PM10      PM2_5      NMVOC NH3 Unit Quantty Type
## 1   11.5790  0.6673455  0.09684276 0.03978916 0.08637947   0   NA      NA   NA
## 2   72.5170  4.1750000  1.06799662 0.43880094 0.54097761   0   NA      NA   NA
## 3  284.5313 36.2251000  3.57381777 1.46835164 2.12260659   0   NA      NA   NA
## 4   27.4530  1.1820151  0.19795276 0.08133158 0.20479968   0   NA      NA   NA
## 5  112.1674 27.8750000  2.25393515 0.92605991 0.83677002   0   NA      NA   NA
## 6   21.2126  1.2580000  0.42127192 0.17308530 0.15824623   0   NA      NA   NA
## 7   23.1560  0.9121350  0.21193984 0.08707837 0.17274401   0   NA      NA   NA
## 8   39.2910 96.0010000 18.70001630 7.68315603 0.29311129   0   NA      NA   NA
## 9   51.9012  0.3656000  0.39069000 0.16052030 0.38718352   0   NA      NA   NA
## 10  23.6490  2.5205763  0.17792041 0.07310102 0.17642180   0   NA      NA   NA
##                                   X___10            Region         Area
## 1   Beogradske elektrane, TO Banovo brdo Beogradski region Grad Beograd
## 2         Beogradske elektrane, TO Cerak Beogradski region Grad Beograd
## 3  Beogradske elektrane, TO Novi Beograd Beogradski region Grad Beograd
## 4    Beogradske elektrane, TO Miljakovac Beogradski region Grad Beograd
## 5         Beogradske elektrane, TO Dunav Beogradski region Grad Beograd
## 6     Beogradske elektrane, TO Medaković Beogradski region Grad Beograd
## 7      Beogradske elektrane, TO Voždovac Beogradski region Grad Beograd
## 8         Beogradske elektrane, TO Zemun Beogradski region Grad Beograd
## 9     Beogradske elektrane, TO Konjarnik Beogradski region Grad Beograd
## 10     Beogradske elektrane, TO Mirijevo Beogradski region Grad Beograd
##                 Mncplty                Locatin                 geometry
## 1      Beograd-Čukarica     Beograd (Čukarica) POINT (453732.9 4958584)
## 2      Beograd-Čukarica     Beograd (Čukarica) POINT (453736.5 4953799)
## 3  Beograd-Novi Beograd Beograd (Novi Beograd) POINT (454145.1 4961030)
## 4      Beograd-Rakovica     Beograd (Rakovica) POINT (456831.2 4955267)
## 5    Beograd-Stari Grad   Beograd (Stari Grad)   POINT (456269 4963981)
## 6      Beograd-Voždovac     Beograd (Voždovac) POINT (459853.8 4958226)
## 7      Beograd-Voždovac     Beograd (Voždovac) POINT (459553.2 4955717)
## 8         Beograd-Zemun        Beograd (Zemun) POINT (451897.5 4965798)
## 9      Beograd-Zvezdara     Beograd (Zvezdara) POINT (453732.9 4958584)
## 10     Beograd-Zvezdara     Beograd (Zvezdara) POINT (463136.6 4959595)

sf.1A1a_ep %<>% st_transform(4326)

Učitavanje linijskih izvora zagađivača i transformacija:

sf.1A3bi_R <- st_read(dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/1A3bi_R.shp")

## Reading layer `1A3bi_R' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\1A3bi_R.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1440 features and 23 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 333615.1 ymin: 4656734 xmax: 657521.3 ymax: 5114098
## projected CRS:  WGS 84 / UTM zone 34N

sf.1A3bi_R

## Simple feature collection with 1440 features and 23 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 333615.1 ymin: 4656734 xmax: 657521.3 ymax: 5114098
## projected CRS:  WGS 84 / UTM zone 34N
## First 10 features:
##    ID is_PGDS PGDS_2015_ Name Kategrj Broj_pt Oznak_d Smer Oznk_pc Oznk_zv
## 1   1       0   1504.233  IIA     IIA     100   10001    O   10001    1301
## 2   2       1   1859.667  IIA     IIA     100   10002    O    1302   10002
## 3   3       0   4787.406  IIA     IIA     100   10003    O   10002     102
## 4   4       1   4787.371  IIA     IIA     100 10004o1    O     102   102.1
## 5   5       0   6670.213  IIA     IIA     100 10004o2    D   102.1   10003
## 6   6       0   6669.752  IIA     IIA     100 10004o3    L   10003   102.1
## 7   7       0   6670.213  IIA     IIA     100 10005o1    O   10003 10003.1
## 8   8       0   6666.552  IIA     IIA     100 10005o2    D 10003.1 10003.2
## 9   9       0   6670.213  IIA     IIA     100 10005o3    L 10003.2 10003.1
## 10 10       0   4273.142  IIA     IIA     100 10005o4    O 10003.2    1103
##    Duzina Stacinz                           Nazv_pc               Nzv_zvr
## 1   4.679   4.679 granica MA<U+00D0>/SRB (Horgoš 2)       Horgoš (put A1)
## 2  10.439  15.118                  Horgoš (Kanjiža)       Backi Vinogradi
## 3     3.5  18.618                   Backi Vinogradi petlja Subotica sever
## 4   2.349  20.967             petlja Subotica sever                <Null>
## 5   9.265  30.232                            <Null>     Subotica (centar)
## 6   9.244  <Null>                 Subotica (centar)                <Null>
## 7   0.983  31.215                 Subotica (centar)                <Null>
## 8   2.394  33.609                            <Null>                <Null>
## 9   2.402  <Null>                            <Null>                <Null>
## 10  2.934  36.543                            <Null>   Subotica (B.Topola)
##    Komentr PGDS_2015       NOx        SO2       PM10      PM2_5     NMVOC
## 1   <Null>  1504.233 0.8676208 0.03817683 0.05566276 0.04222431 0.3049724
## 2   <Null>  1377.000 0.7942347 0.03494771 0.05095463 0.03865285 0.2791768
## 3   <Null>  4787.406 2.7613098 0.12150244 0.17715357 0.13438407 0.9706120
## 4   <Null>  6673.000 3.8488946 0.16935807 0.24692826 0.18731333 1.3529026
## 5   <Null>  6670.213 3.8472872 0.16928734 0.24682514 0.18723510 1.3523376
## 6   <Null>  6669.752 3.8470213 0.16927564 0.24680808 0.18722216 1.3522442
## 7   <Null>  6670.213 3.8472872 0.16928734 0.24682514 0.18723510 1.3523376
## 8   <Null>  6666.552 3.8451755 0.16919442 0.24668966 0.18713233 1.3515953
## 9   <Null>  6670.213 3.8472872 0.16928734 0.24682514 0.18723510 1.3523376
## 10  <Null>  4273.142 2.4646895 0.10845063 0.15812371 0.11994851 0.8663487
##           NH3     PGDS                       geometry
## 1  0.03728672 1504.233 LINESTRING (418315.4 511096...
## 2  0.03413289 1377.000 LINESTRING (410263.6 510721...
## 3  0.11866955 4787.406 LINESTRING (406835.8 510659...
## 4  0.16540940 6673.000 LINESTRING (404627.4 510609...
## 5  0.16534032 6670.213 LINESTRING (398266.3 510462...
## 6  0.16532890 6669.752 LINESTRING (404627.1 510609...
## 7  0.16534032 6670.213 LINESTRING (397392 5104191,...
## 8  0.16524957 6666.552 LINESTRING (397840.2 510196...
## 9  0.16534032 6670.213 LINESTRING (397394.8 510419...
## 10 0.10592205 4273.142 LINESTRING (396709.4 509936...

sf.1A3bi_R %<>% st_transform(4326)

Učitavanje površinskih izvora zagađivača - pomoćni podaci i transformacija:

CLC_12_18 <- st_read(dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/CLC18_RS.shp")

## Reading layer `CLC18_RS' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\CLC18_RS.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34794 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 331043.4 ymin: 4674468 xmax: 663570 ymax: 5117030
## projected CRS:  ETRS89 / UTM zone 34N

CLC_12_18

## Simple feature collection with 34794 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 331043.4 ymin: 4674468 xmax: 663570 ymax: 5117030
## projected CRS:  ETRS89 / UTM zone 34N
## First 10 features:
##       ID CODE_18  Area_Ha Remark SHAPE_Leng SHAPE_Area
## 1   RS_1     111 57.65009   <NA>   6456.549   576500.9
## 2   RS_2     111 37.62057   <NA>   5540.471   376205.7
## 3   RS_3     111 35.45603   <NA>   3725.036   354560.3
## 4   RS_4     111 25.14260   <NA>   4178.823   251426.0
## 5   RS_5     112 42.43190   <NA>   4099.389   424319.0
## 6   RS_6     112 67.68868   <NA>   4377.618   676886.8
## 7   RS_7     112 28.49615   <NA>   2736.777   284961.5
## 8   RS_8     112 25.46822   <NA>   2175.572   254682.2
## 9   RS_9     112 39.81418   <NA>   3961.963   398141.8
## 10 RS_10     112 34.70737   <NA>   2637.301   347073.7
##                          geometry
## 1  POLYGON ((574615.8 4796652,...
## 2  POLYGON ((572794.7 4796959,...
## 3  POLYGON ((457109.9 4962932,...
## 4  POLYGON ((397083.5 5106354,...
## 5  POLYGON ((555111.5 4676318,...
## 6  POLYGON ((555062.8 4679080,...
## 7  POLYGON ((560340.8 4680704,...
## 8  POLYGON ((553657.7 4681341,...
## 9  POLYGON ((559394.7 4682410,...
## 10 POLYGON ((553051 4682493, 5...

CLC_12_18 %<>% st_transform(4326)

---

Učitavanje grida i transformacija u koordinatni sistem u projekciji (UTM projekcija, 34N zona / WGS84 elipsoid):

gridWGS84 <- st_read(dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/Grid_Serbia_0.05deg.gpkg")

## Reading layer `Grid_Serbia_005' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\Grid_Serbia_0.05deg.gpkg' using driver `GPKG'
## Simple feature collection with 3769 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 18.8 ymin: 42.2 xmax: 23 ymax: 46.2
## geographic CRS: WGS 84

gridWGS84

## Simple feature collection with 3769 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 18.8 ymin: 42.2 xmax: 23 ymax: 46.2
## geographic CRS: WGS 84
## First 10 features:
##     left   top right bottom  id value                           geom
## 1  18.80 45.95 18.85  45.90   6     1 POLYGON ((18.8 45.95, 18.85...
## 2  18.80 45.90 18.85  45.85   7     1 POLYGON ((18.8 45.9, 18.85 ...
## 3  18.80 45.85 18.85  45.80   8     1 POLYGON ((18.8 45.85, 18.85...
## 4  18.85 45.95 18.90  45.90  96     1 POLYGON ((18.85 45.95, 18.9...
## 5  18.85 45.90 18.90  45.85  97     1 POLYGON ((18.85 45.9, 18.9 ...
## 6  18.85 45.85 18.90  45.80  98     1 POLYGON ((18.85 45.85, 18.9...
## 7  18.85 45.80 18.90  45.75  99     1 POLYGON ((18.85 45.8, 18.9 ...
## 8  18.85 45.75 18.90  45.70 100     1 POLYGON ((18.85 45.75, 18.9...
## 9  18.85 45.65 18.90  45.60 102     1 POLYGON ((18.85 45.65, 18.9...
## 10 18.85 45.60 18.90  45.55 103     1 POLYGON ((18.85 45.6, 18.9 ...

Presek i sumiranje vrednosti polutanata po kategorijama po ćelijama grida

---

1A1a - Public electricity production (tačkasti izvor)

---

vars <- c("NOx","SO2", "PM10", "PM2_5", "NMVOC", "NH3")
sf.1A1a_ep %<>% dplyr::select(vars)
names(sf.1A1a_ep)

## [1] "NOx"      "SO2"      "PM10"     "PM2_5"    "NMVOC"    "NH3"      "geometry"

Prostorizacija - sumiranje vrednosti tačkastih zagađivača po ćelijama grida:

p.1A1a_ep <- gridWGS84 %>%
  sf::st_join(sf.1A1a_ep) %>%
  dplyr::group_by(id) %>%
  dplyr::summarize(NOx = sum(NOx, na.rm = TRUE),
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2_5 = sum(PM2_5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  dplyr::mutate(id = as.numeric(id))

Kontrola:

sum.sf.1A1a_ep <- sf.1A1a_ep %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)

sum.p.1A1a_ep <- p.1A1a_ep %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)


data.frame(sum = c("prostorizovano", "ukupno", "razlika"), rbind(sum.p.1A1a_ep, sum.sf.1A1a_ep, data.frame(sum.p.1A1a_ep -  sum.sf.1A1a_ep))) %>%
 datatable(caption = "Tabela 1: Sumarne vrednosti razlika nakon prostorizacije") 

Vizuelizacija:

map.1A1a_ep <- p.1A1a_ep
map.1A1a_ep$Spatialised <- 0
map.1A1a_ep$Spatialised[map.1A1a_ep$NOx != 0 | map.1A1a_ep$SO2 != 0 | map.1A1a_ep$PM10 != 0 | map.1A1a_ep$PM2_5 != 0 | map.1A1a_ep$NMVOC != 0 | map.1A1a_ep$NH3 != 0] <- 1 # Zadatak: resiti sa mutate i case_when

mapview(map.1A1a_ep, zcol = "Spatialised", layer.name = "Spatialised 1A1a_ep")

1A3bi_R - Road Transport: Passengers cars - Rural transport (linijski izvor)

---

Mreža puteva po kategorijama puteva:

ggplot()+
  geom_sf(data = sf.1A3bi_R, aes(colour = Name))+
  scale_colour_manual(name = "Kategorija: ", values = c("red", "orange"))+
  labs(title = "Mreža puteva", 
       subtitle = "Klasifikacija po kategorijama puteva",
       caption = "Republika Srbija")+
  theme_bw()+
  theme(legend.position = "bottom")

Mreža puteva po vrednosti protoka saobraćaja:

ggplot()+
  geom_sf(data = sf.1A3bi_R, aes(size = PGDS_2015/10000), colour = "black")+
  scale_size_identity()+
  labs(title = "Mreža puteva", 
       subtitle = "Klasifikacija po vrednosti protoka saobraćaja",
       caption = "Republika Srbija")+
  theme_bw()+
  theme(legend.position = "bottom")

sf.1A3bi_R %<>% dplyr::select(ID, vars, PGDS_2015) 

Sumiranje ukupnih vrednosti po polutantima:

sum.sf.1A3bi_R <- sf.1A3bi_R %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)
sum.sf.1A3bi_R %>% datatable()

Presek sa gridom:

sf.1A3bi_R %<>% dplyr::mutate(Length = st_length(.),
                              PGDSL = PGDS_2015 * Length)

sf.1A3bi_R.int <- st_intersection(sf.1A3bi_R, gridWGS84) %>%
  dplyr::select(., PGDS_2015, PGDSL, Length) %>%
  mutate(Length_int = st_length(.), 
         PGDS_int = (PGDS_2015/Length)*Length_int) %>%
  dplyr::select(., Length, PGDS_int)
  
sum_PGDS <- sum(sf.1A3bi_R.int$PGDS_int)

sf.1A3bi_R.int <- sf.1A3bi_R.int %>%
  dplyr::mutate(NOx = ((sum.sf.1A3bi_R$NOx/sum_PGDS)*PGDS_int),
         SO2 = ((sum.sf.1A3bi_R$SO2/sum_PGDS)*PGDS_int),
         PM10 = ((sum.sf.1A3bi_R$PM10/sum_PGDS)*PGDS_int),
         PM2_5 = ((sum.sf.1A3bi_R$PM2_5/sum_PGDS)*PGDS_int),
         NMVOC = ((sum.sf.1A3bi_R$NMVOC/sum_PGDS)*PGDS_int),
         NH3 = ((sum.sf.1A3bi_R$NH3/sum_PGDS)*PGDS_int)) %>% 
  dplyr::select(vars)

# Kontrola

sum.sf.1A3bi_R.int <- sf.1A3bi_R.int %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)

sum.sf.1A3bi_R

##       NOx    SO2   PM10  PM2_5   NMVOC    NH3
## 1 3099.65 136.39 198.86 150.85 1089.54 133.21

sum.sf.1A3bi_R.int

##       NOx    SO2   PM10  PM2_5   NMVOC    NH3
## 1 3099.65 136.39 198.86 150.85 1089.54 133.21

Prostorizacija - sumiranje vrednosti linijskih zagađivača po ćelijama grida:

p.1A3bi_R <- gridWGS84 %>%
  st_join(sf.1A3bi_R.int, join = st_contains) %>% 
  mutate(id = as.numeric(id)) %>%
  group_by(id) %>%
  summarize(NOx = sum(NOx, na.rm = TRUE),
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2_5 = sum(PM2_5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  mutate(id = as.numeric(id))
p.1A3bi_R[,vars] <- drop_units(p.1A3bi_R[,vars])

Kontrola:

sum.sf.1A3bi_R <- sf.1A3bi_R %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)

sum.p.1A3bi_R <- p.1A3bi_R %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)


data.frame(sum = c("prostorizovano", "ukupno", "razlika"), rbind(sum.p.1A3bi_R, sum.sf.1A3bi_R, data.frame(sum.p.1A3bi_R -  sum.sf.1A3bi_R))) %>%
 datatable(caption = "Tabela 2: Sumarne vrednosti razlika nakon prostorizacije") 

Vizuelizacija:

map.1A3bi_R <- p.1A3bi_R
map.1A3bi_R$Spatialised <- 0
#map.1A3bi_R[,vars] <- units::drop_units(map.1A3bi_R[,vars])
map.1A3bi_R$Spatialised[map.1A3bi_R$NOx != 0 | map.1A3bi_R$SO2 != 0 | map.1A3bi_R$PM10 != 0 | map.1A3bi_R$PM2_5 != 0 | map.1A3bi_R$NMVOC != 0 | map.1A3bi_R$NH3 != 0] <- 1 

mapview(map.1A3bi_R, zcol = "Spatialised", layer.name = "Spatialised 1A3bi_R")

3Da1 - Inorganic N-fertilizers (površinski izvor)

---

Corine Land Cover karta (zemljišni pokrivač):

ggplot()+
  geom_sf(data = CLC_12_18, aes(fill = CODE_18))+
  #scale_colour_manual(name = "Kategorija: ", values = c("red", "orange"))+
  labs(title = "Corine Land Cover karta", 
       subtitle = "Klasifikacija po kategorijama zemljišnog pokrivača",
       caption = "Republika Srbija")+
  theme_bw()+
  theme(legend.position = "right")

Ukupne vrednosti po polutantima:

sum.sf.3Da1 <- data.frame(NOx = 8.33, SO2 = 0, PM10 = 0, PM2_5 = 0, NMVOC = 0, NH3 = 18028.55)
sum.sf.3Da1

##    NOx SO2 PM10 PM2_5 NMVOC      NH3
## 1 8.33   0    0     0     0 18028.55

Izdvajanje poljoprivrednog zemljišta (CLC classes) i presek sa gridom:

sf_clc18_polj <- subset(CLC_12_18, CODE_18 == "211" | CODE_18 == "221" | CODE_18 == "222" | CODE_18 == "242" | CODE_18 == "243") # CLC agricultural areas
  
sf_clc18_polj[,vars] <- NA

Presek sa gridom:

sf_clc18_polj.int <- st_intersection(sf_clc18_polj, gridWGS84) %>%
  dplyr::mutate(Area = st_area(.)) # racunanje povrsina poligona

sf_clc18_polj.int %<>% dplyr::select(., vars, Area) 

Interpolacija vrednosti po polutantima srazmerno površini entiteta:

sum_Area <- sum(sf_clc18_polj.int$Area)

sf.3Da1 <- sf_clc18_polj.int %>%
  mutate(NOx = ((sum.sf.3Da1$NOx/sum_Area)*Area),
         SO2 = ((sum.sf.3Da1$SO2/sum_Area)*Area),
         PM10 = ((sum.sf.3Da1$PM10/sum_Area)*Area),
         PM2_5 = ((sum.sf.3Da1$PM2_5/sum_Area)*Area),
         NMVOC = ((sum.sf.3Da1$NMVOC/sum_Area)*Area),
         NH3 = ((sum.sf.3Da1$NH3/sum_Area)*Area))
sf.3Da1 %<>% dplyr::select(vars)

Prostorizacija - sumiranje vrednosti površinskih zagađivača po ćelijama grida:

p.3Da1 <- gridWGS84 %>%
  st_join(sf.3Da1, join = st_contains) %>% 
  group_by(id) %>%
  summarize(NOx = sum(NOx, na.rm = TRUE),
            SO2 = sum(SO2, na.rm = TRUE),
            PM10 = sum(PM10, na.rm = TRUE),
            PM2_5 = sum(PM2_5, na.rm = TRUE),
            NMVOC = sum(NMVOC, na.rm = TRUE),
            NH3 = sum(NH3, na.rm = TRUE)) %>% 
  mutate(id = as.numeric(id))

p.3Da1[, vars] <- drop_units(p.3Da1[, vars])

Kontrola:

sum.p.3Da1 <- p.3Da1 %>%
  sf::st_drop_geometry() %>%
  dplyr::select(., vars) %>% 
  apply(., 2, sum) %>% 
  t(.) %>% 
  as.data.frame() %>%
  dplyr::mutate_if(is.numeric, round, 2)


data.frame(sum = c("prostorizovano", "ukupno", "razlika"), rbind(sum.p.3Da1, sum.sf.3Da1, data.frame(sum.p.3Da1 -  sum.sf.3Da1))) %>%
 datatable(caption = "Tabela 2: Sumarne vrednosti razlika nakon prostorizacije") 

Vizuelizacija:

map.3Da1 <- p.3Da1
map.3Da1$Spatialised <- 0
#map.3Da1[,vars] <- units::drop_units(map.3Da1[,vars])
map.3Da1$Spatialised[map.3Da1$NOx != 0 | map.3Da1$SO2 != 0 | map.3Da1$PM10 != 0 | map.3Da1$PM2_5 != 0 | map.1A3bi_R$NMVOC != 0 | map.3Da1$NH3 != 0] <- 1 

mapview(map.3Da1, zcol = "Spatialised", layer.name = "Spatialised 1A3bi_R")

Kreiranje jedinstvenog grida

---

data.spat.list <- list(p.1A1a_ep, p.1A3bi_R, p.3Da1)                                                   

sf_data <- data.spat.list[[1]]
for(i in 2:length(data.spat.list)){                                             
 sf_data <- st_join(sf_data, data.spat.list[[i]], join = st_equals) %>% 
   group_by(id.x) %>%
   summarize(NOx = sum(NOx.x, NOx.y),
             SO2 = sum(SO2.x, SO2.y),
             PM10 = sum(PM10.x, PM10.y),
             PM2_5 = sum(PM2_5.x, PM2_5.y),
             NMVOC = sum(NMVOC.x, NMVOC.y),
             NH3 = sum(NH3.x + NH3.y)) %>%
   mutate(id = id.x) %>%
   select(id, NOx, SO2, PM10, PM2_5, NMVOC, NH3)
 print(paste("NOx:",sum(sf_data$NOx))) 
 print(paste("SO2:",sum(sf_data$SO2)))
 print(paste("PM10:",sum(sf_data$PM10)))
 print(paste("PM2_5:",sum(sf_data$PM2_5)))
 print(paste("NMVOC:",sum(sf_data$NMVOC)))
 print(paste("NH3:",sum(sf_data$NH3)))
}

## [1] "NOx: 4339.7786"
## [1] "SO2: 1910.84306396663"
## [1] "PM10: 323.042992264955"
## [1] "PM2_5: 201.872271321156"
## [1] "NMVOC: 1098.79137281681"
## [1] "NH3: 133.21"
## [1] "NOx: 4348.1086"
## [1] "SO2: 1910.84306396663"
## [1] "PM10: 323.042992264955"
## [1] "PM2_5: 201.872271321156"
## [1] "NMVOC: 1098.79137281681"
## [1] "NH3: 18161.76"

sf_data

## Simple feature collection with 3769 features and 7 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 18.8 ymin: 42.2 xmax: 23 ymax: 46.2
## geographic CRS: WGS 84
## # A tibble: 3,769 x 8
##       id     NOx     SO2    PM10   PM2_5  NMVOC    NH3                      geom
##    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>    <POLYGON [arc_degree]>
##  1     6 0.      0       0       0       0      0      ((18.8 45.95, 18.85 45.9~
##  2     7 4.25e-5 0       0       0       0      0.0920 ((18.8 45.9, 18.85 45.9,~
##  3     8 1.84e-4 0       0       0       0      0.398  ((18.8 45.85, 18.85 45.8~
##  4    96 1.23e-4 0       0       0       0      0.266  ((18.85 45.95, 18.9 45.9~
##  5    97 2.16e-1 0.00941 0.0137  0.0104  0.0751 4.31   ((18.85 45.9, 18.9 45.9,~
##  6    98 1.08e-1 0.00470 0.00686 0.00520 0.0376 1.70   ((18.85 45.85, 18.9 45.8~
##  7    99 0.      0       0       0       0      0      ((18.85 45.8, 18.9 45.8,~
##  8   100 0.      0       0       0       0      0      ((18.85 45.75, 18.9 45.7~
##  9   102 0.      0       0       0       0      0      ((18.85 45.65, 18.9 45.6~
## 10   103 0.      0       0       0       0      0      ((18.85 45.6, 18.9 45.6,~
## # ... with 3,759 more rows

Vizuelizacija

---

classes.NOx <- classIntervals(sf_data$NOx, n = 12, style = "fisher")
classes.SO2 <- classIntervals(sf_data$SO2, n = 12, style = "fisher")
classes.PM10 <- classIntervals(sf_data$PM10, n = 12, style = "fisher")
classes.PM2_5 <- classIntervals(sf_data$PM2_5, n = 12, style = "fisher")
classes.NMVOC <- classIntervals(sf_data$NMVOC, n = 12, style = "fisher")
classes.NH3 <- classIntervals(sf_data$NH3, n = 12, style = "fisher")

sf_data <- sf_data %>%
  mutate(percent_class_NOx = cut(NOx, classes.NOx$brks, include.lowest = T,dig.lab=5),
         percent_class_SO2 = cut(SO2, classes.SO2$brks, include.lowest = T,dig.lab=7),
         percent_class_PM10 = cut(PM10, classes.PM10$brks, include.lowest = T,dig.lab=5),
         percent_class_PM2_5 = cut(PM2_5, classes.PM2_5$brks, include.lowest = T,dig.lab=5),
         percent_class_NMVOC = cut(NMVOC, classes.NMVOC$brks, include.lowest = T,dig.lab=5),
         percent_class_NH3 = cut(NH3, classes.NH3$brks, include.lowest = T,dig.lab=5)
         )

pal1 <- viridisLite::viridis(12, direction = -1)

a<-ggplot() +
  geom_sf(data = sf_data,
          aes(fill = percent_class_NOx)) +
  scale_fill_manual(values = pal1,
                    name = "NOx [t]") +
  
  labs(x = NULL, y = NULL,
       title = "Pollutant inventory spatialization - NOx",
       subtitle = "Spatial resolution 0.05°x0.05°",
       caption = "UBFCE (2020)") +
  theme(line = element_blank(),
        #axis.text = element_blank(),
        axis.title = element_blank(),
        #legend.position = "None", ###################### legend
        panel.background = element_blank()) +
  coord_sf(datum = sf::st_crs(4326))

b<-ggplot() +
  geom_sf(data = sf_data,
          aes(fill = percent_class_SO2)) +
  scale_fill_manual(values = pal1,
                    name = "SO2 [t]") +
  labs(x = NULL, y = NULL,
       title = "Pollutant inventory spatialization - SO2",
       subtitle = "Spatial resolution 0.05°x0.05°",
       caption = "UBFCE (2020)") +
  theme(line = element_blank(),
        #axis.text = element_blank(),
        axis.title = element_blank(),
        #legend.position = "None", ###################### legend
        panel.background = element_blank()) +
  coord_sf(datum = sf::st_crs(4326))

c<-ggplot() +
  geom_sf(data = sf_data,
          aes(fill = percent_class_PM10)) +
  scale_fill_manual(values = pal1,
                    name = "PM10 [t]") +
  labs(x = NULL, y = NULL,
       title = "Pollutant inventory spatialization - PM10",
       subtitle = "Spatial resolution 0.05°x0.05°",
       caption = "UBFCE (2020)") +
  theme(line = element_blank(),
        #axis.text = element_blank(),
        axis.title = element_blank(),
        #legend.position = "None", ###################### legend
        panel.background = element_blank()) +
  coord_sf(datum = sf::st_crs(4326))

d<-ggplot() +
  geom_sf(data = sf_data,
          aes(fill = percent_class_PM2_5)) +
  scale_fill_manual(values = pal1,
                    name = "PM2.5 [t]") +
  labs(x = NULL, y = NULL,
       title = "Pollutant inventory spatialization - PM2.5",
       subtitle = "Spatial resolution 0.05°x0.05°",
       caption = "UBFCE (2020)") +
  theme(line = element_blank(),
        #axis.text = element_blank(),
        axis.title = element_blank(),
        #legend.position = "None", ###################### legend
        panel.background = element_blank()) +
  coord_sf(datum = sf::st_crs(4326))

e<-ggplot() +
  geom_sf(data = sf_data,
          aes(fill = percent_class_NMVOC)) +
  scale_fill_manual(values = pal1,
                    name = "NMVOC [t]") +
  labs(x = NULL, y = NULL,
       title = "Pollutant inventory spatialization - NMVOC",
       subtitle = "Spatial resolution 0.05°x0.05°",
       caption = "UBFCE (2020)") +
  theme(line = element_blank(),
        #axis.text = element_blank(),
        axis.title = element_blank(),
        #legend.position = "None", ###################### legend
        panel.background = element_blank()) +
  coord_sf(datum = sf::st_crs(4326))

f<-ggplot() +
  geom_sf(data = sf_data,
          aes(fill = percent_class_NH3)) +
  scale_fill_manual(values = pal1,
                    name = "NH3 [t]") +
  labs(x = NULL, y = NULL,
       title = "Pollutant inventory spatialization - NH3",
       subtitle = "Spatial resolution 0.05°x0.05°",
       caption = "UBFCE (2020)") +
  theme(line = element_blank(),
        #axis.text = element_blank(),
        axis.title = element_blank(),
        #legend.position = "None", ###################### legend
        panel.background = element_blank()) +
  coord_sf(datum = sf::st_crs(4326))

grid.arrange(a, b, c, d, e, f, ncol = 2, nrow = 3)

Čuvanje podataka u geopackage formatu za razmenu prostornih podataka

st_write(sf_data, dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/grid_final.gpkg", layer = "Final_grid")

st_write(p.1A1a_ep, dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/grid_final.gpkg", layer = "1A1a_ep - Public electricity production")

st_write(p.1A3bi_R, dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/grid_final.gpkg", layer = "1A3bi_R - Road Transport: Passengers cars - Rural transport")

st_write(p.3Da1, dsn = "D:/R_projects/Nauka_R/Slides/R/Spatial data vezbe/Data/grid_final.gpkg", layer = "3Da1 - Inorganic N-fertilizers")